#ifndef _NEWTON_INTERPOLATION_H
#define _NEWTON_INTERPOLATION_H

#include <iostream>
#include <math.h>
#include <vector>
#include "Polynomial.h"

using namespace std;

int fact(int n)
{
    if(n > 1)
        return n * fact(n-1);
    else
        return 1;
}

class Newton 
{
    protected:
    std::vector<double> x;
    std::vector<double> f;
    public:
    Polynomial poly;
    int n = 0;
    Newton(){}
    Newton(std::vector<double> _x,std::vector<double> _f)
    {
        x = _x;
        f = _f;
        n = x.size() - 1;
    }
    ~Newton()
    {
        x.clear();
        f.clear();
    }
    
    Polynomial polation()
    {
        Polynomial p(f[0]);
        Polynomial monos(1);    //单项式连乘
        std::vector<double> old(f);
        for (int i = 0; i < n; i++)
        {
            std::vector<double> temp(n-i);
            for (int j = 0; j < n-i; j++)
            {
                temp[j] = (old[j+1] - old[j]) / (x[j+1+i] - x[j]);
            }
            Polynomial mono(-1*x[i],1);     //单项式
            monos = monos * mono;
            p = p + (monos * temp[0]);
            old = temp;

            /*
            for (int k = 0; k < temp.size(); k++)
            {
                cout << temp[k] << "\t";
            }
            cout << endl;*/


        }
        poly = p;
        return p;
    }

};

//Hermite插值能覆盖Newton插值
class Hermite : public Newton
{
    public:
    Hermite(std::vector<double> _x,std::vector<double> _f)
    {
        x = _x;
        f = _f;
        n = x.size() - 1;
    }
    ~Hermite()
    {
        x.clear();
        f.clear();
    }
    
    Polynomial polation()
    {
        Polynomial p(f[0]);
        Polynomial monos(1);    //单项式连乘
        std::vector<double> old(f);
        for (int i = 1; i <= n; i++)
        {
            if (x[i] == x[i-1])
            {
                old[i] = old[i-1];
            }
        }
        for (int i = 0; i < n; i++)
        {
            std::vector<double> temp(n-i);
            int cnt = 0;
            for (int j = 0; j < n-i; j++)
            {
                if (x[j+i+1] == x[j])
                {
                    temp[j] = f[j+i+1-cnt] / fact(i+1);
                    cnt++;
                }
                else
                {
                    cnt = 0;
                    temp[j] = (old[j+1] - old[j]) / (x[j+1+i] - x[j]);
                }
            }
            Polynomial mono(-1*x[i],1);     //单项式
            monos = monos * mono;
            p = p + (monos * temp[0]);
            old = temp;
            /*for (int k = 0; k < temp.size(); k++)
            {
                cout << temp[k] << "\t";
            }
            cout << endl;*/
        }
        poly = p;
        return p;
    }

};



#else
    ///do nothing
#endif